import matplotlib.pyplot as plt
import numpy as np
import glob

plt.rc('axes', titlesize=14)     # fontsize of the axes title
plt.rc('axes', labelsize=14)    # fontsize of the x and y labels
plt.rc('xtick', labelsize=14)    # fontsize of the tick labels
plt.rc('ytick', labelsize=14)    # fontsize of the tick labels

fl=glob.glob('capture-results/instability/*.osc')
print(fl)

for f in fl:
    a,e,i=np.loadtxt(f,unpack=True,usecols=(1,2,3))
    ind=np.argwhere(i<15)
    plt.scatter(a[ind],e[ind],c='#888888',s=5)

Athor=[2.38,  0.1,  0.15]
AthorReg=(2.3, 2.48)

x=np.linspace(2,2.5,100)
y=1-1.8/x
plt.plot(x,y,c='k')
#y=1-2.0/x
#plt.plot(x,y,c='#BBBBBB')


x=np.linspace(AthorReg[0],AthorReg[1])
yUp=1-2.0/x
yDo=np.zeros(len(x))
plt.fill_between(x, yUp, yDo, color='#BBBBBB',
                 alpha=0.5)

plt.scatter(Athor[0],Athor[1],s=30,c='k',zorder=10)

plt.ylim(0,1)
plt.xlim(2.0, 2.5)
plt.text(2.364, 0.055, "Athor")
plt.xlabel('Semimajor axis, $a$ (au)')
plt.ylabel('Eccentricity, $e$')
plt.tight_layout()
plt.savefig('Fig3.pdf', dpi=300)
